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Abstract 

We introduce a new methodology for adding localized, space-time smooth, artificial viscosity to 
nonlinear systems of conservation laws which propagate shock waves, rarefactions, and contact 
discontinuities, which we call the C-method. We shall focus our attention on the compressible Euler 
equations in one space dimension. The novel feature of our approach involves the coupling of a linear 
scalar reaction-diffusion equation to our system of conservation laws, whose solution C(.t, t) is the 
coefficient to an additional (and artificial) term added to the flux, which determines the location, 
localization, and strength of the artificial viscosity. Near shock discontinuities, C{x,t) is large and 
localized, and transitions smoothly in space-time to zero away from discontinuities. Our approach 
is a provably convergent, spacetime-regularized variant of the original idea of Richtmeyer and Von 
Neumann, and is provided at the level of the PDE, thus allowing a host of numerical discretization 
schemes to be employed. 

We demonstrate the effectiveness of the C-method with three different numerical implementations 
and apply these to a collection of classical problems: the Sod shock-tube, the Osher-Shu shock-tube, 
the Woodward- Colella blast wave and the Lcblanc shock-tube. First, we use a classical continuous 
finite-element implementation using second-order discretization in both space and time , FEM-C. 
Second, we use a simplified WENO scheme within our C-method framework, WENO-C. Third, we 
use WENO with the Lax-Friedrichs flux together with the C-equation, and call this WENO-LF-C. 
All three schemes yield higher-order discretization strategies, which provide sharp shock resolution 
with minimal overshoot and noise, and compare well with higher-order WENO schemes that employ 
approximate Riemann solvers, outperforming them for the difficult Leblanc shock tube experiment. 
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1. Introduction 



1.1. Smoothing conservation laws 

The initial-value problem for a general nonlinear system of conservation laws can be written as 
an evolution equation, 

dtU{x, t) + div F{U{x, t)) = with U\t=o = Uq , (1) 

for an m-vector U defined on (Z?+l)-dimensional space-time. Such partial differential equations 
(PDE) are both ubiquitous and fundamental in science and engineering, and include the compressible 
Euler equations of gas dynamics, the magneto-hydrodynamic (MHD) equations modeling ionized 
plasma, the elasticity equations of solid mechanics, and numerous related physical systems which 
possess complicated nonlinear wave interactions. 

It is well known that solutions of ([!]) can develop finite-time shocks, even when the initial data 
is smooth, in which case, discontinuities of U are propagated according to the so-called Rankine- 
Hugoniot conditions (see Section [2 . 1 1 below) . It is important to develop stable and robust numerical 
algorithms which can approximate weak shock-wave solutions. Even in one-space dimension, non- 
linear wave interaction such as two shock waves colliding, is a difficult problem when considering 
accuracy, stability and monotonicity. The challenge is maintaining higher-order accuracy away from 
the shock while approximating the discontinuity in an order-Ao; smooth transition region where Ax 
denotes the spatial grid size. 

As we describe below, a variety of clever discretization schemes have been developed and em- 
ployed, particularly in one-space dimension, to approximate discontinuous solution profiles in an 
essentially non-oscillatory (ENO) fashion. These include, but are not limited to, total variation 
diminishing (TVD) schemes, flux-corrected transport (FCT) schemes, weighted essentially non- 
oscillatory (WENO) schemes, discontinuous Galerkin methods, artificial diffusion methods, exact 
and approximate Riemann solvers, and a host of variants and combinations of these techniques. 

We develop a robust parabolic-type regularization of ([T]), which we refer to as the C-method, 
which couples a modified set of m equations for U with an additional linear scalar reaction-diffusion 
equation for a new scalar field C{x, t). Thus, instead of ([T]) we consider a system of m+l equations, 
which use the solution C{x^t) as a coefficient in a carefully chosen modification of the flux. As 
we describe in detail below, the solution C{x,t) is highly localized in regions of discontinuity, and 
transitions smoothly (in both x and t) to zero in regions wherein the solution is smooth. Further, 
as Ax — > 0, we recover the original hyperbolic nonlinear system of conservation laws Q. 

1.2. Numerical discretization 

In the case of 1-D gas dynamics, the construction of non-oscillatory, higher-order, numerical 
algorithms such as ENO by Harten, Engquist, Osher & Chakravarthy [Ij and Shu & Osher [2], 
[3]; WENO by Liu, Osher, & Chan g] and Jiang & Shu 0; MUSCL by Van Leer [6], Colella [7], 
and Huynh or PPM by Colella & Woodward (9l requires carefully chosen reconstruction and 
numerical flux. 

Such numerical methods evolve cell-averaged quantities; to calculate an accurate approximation 
of the flux at cell-interfaces, these schemes reconstruct fcth-order [k > 2) polynomial approximations 
of the solution (and hence the flux) from the computed cell- averages, and thus provide fcth-order 
accuracy away from discontinuities. See, for example, the convergence plots of Greenough & Rider 
[TU] and Liska & Wendroff 11]. Given a polynomial representation of the solution, a strategy 
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is chosen to compute the most accurate cell-interface flux, and this is achieved by a variety of 
algorithms. Centered numerical fluxes, such as Lax-Friedrichs, add dissipation as a mechanism 
to preserve stability and monotonicity. On the other hand, characteristic-type upwinding based 
upon exact (Godunov) or approximate (Roe, Osher, HLL, HLLC) Riemann solvers, which preserve 
monotonicity without adding too much dissipation, tend to be rather complex and PDE-specific; 
moreover, for strong shocks, other techniques may be required to dampen post-shock oscillations 
or to yield entropy-satisfying approximations (see Quirk [H]). Again, we refer the reader to the 
papers [10] , [11] or Colella & Woodward ^3] for a thorough overview, as well as a comparison of the 
effectiveness of a variety of competitive schemes. 

Majda & Osher [14| have shown that any numerical scheme is at best, first-order accurate in 
the presence of shocks or discontinuities. The use of higher-order numerical schemes is, neverthe- 
less, imperative for the elimination of error-terms in the Taylor expansion (in mesh-size) and the 
subsequent limiting of truncation error. Moreover, higher-order schemes tend to be less dissipative 
than there lower-order counterparts, as discussed by Greenough & Rider [lOj : therein, a comparison 
between a 2nd-order PLMDE scheme and a 5th-order WENO scheme demonstrates the improved 
resolution of intricate fine structure afforded by 5th-order WENO, while simultaneously providing 
far less clipping of local extrema than PLMDE. 

In multi-D, similar tools are required to obtain non-oscillatory numerical schemes, but the multi- 
dimensional analogues to those described above are generally limited by mesh considerations. For 
structured grids (such as products of uniform 1-D grids), dimensional splitting is commonly used, 
decomposing the problem into a sequence of 1-D problems. This technique is quite successful, but 
stringent mesh requirements prohibits its use on complex domains. Moreover, applications to PDE 
outside of variants of the Euler equations may be somewhat limited. For further discussion of the 
limitations of dimensional splitting, we refer the reader to Crandall & Majda [H], and Jiang & Tad- 
mor [16j . For unstructured grids, dimensional splitting is not available and alternative approaches 
must be employed, necessitated by the lack of multi-D Riemann solvers. WENO schemes on un- 
structured triangular grids have been developed in Hu & Shu |17| . but using simplified methods, 
which employ reduced characteristic decompositions, can lead to a loss of monotonicity and stability. 

Algorithms that explicitly introduce diffusion provide a simple way to stabilize higher-order nu- 
merical schemes and subsequently remove non-physical oscillations near shocks. In the mathematical 
analysis of conservation laws (and in the truncation error of certain discretization schemes), the sim- 
plest parabolic-regularization is by the addition of a uniform linear viscosity. Choosing a constant 
/3 > 0, which depends upon mesh-size Ax and sometimes velocity or wave-speed, and adding 

f3iA^)dlUix,t) (2) 

to the right hand side of ([l]) provides a uniformly parabolic regularization of the hyperbolic con- 
servation laws, and its discrete implementation smears sharp discontinuities across 0(Ax)-regions 
and thus adds stabilization, but unfortunately, at the cost of accuracy. With the addition of uni- 
form linear viscosity, shocks and discontinuities are captured in a non-oscillatory fashion, but the 
transition region form left to right state, which approximates the discontinuity, tends to grow over 
time. Moreover, since viscosity is applied uniformly over the entire domain I , the benefits of a 
higher-order scheme (away from the discontinuity) may be lost, and the accuracy may often reduce 
to merely first-order. In practical implantation in a numerical scheme, the use of viscosity should 
be localized in regions of shock (so as to stabilize the scheme), limited at contact discontinuities (to 
avoid over-smearing the sharp transition), and very small in smooth regions away from discontinu- 
ities. Achieving these requirements allows higher-order approximation of smooth flow and sharp. 
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non-oscillatory, resolution of shocks and discontinuities. Naturally, this necessitates that the amount 
of added viscosity be a function of the solution. 

The pioneering papers of Richtmyer |T8] , Von Neumann & Richtmyer [19] , Lax & Wendroff [20] , 
and Lapidus [31] suggest the introduction of nonlinear artificial viscosity to equations ([T]) in a form 
similar to the following expression: 

l3{^xfdA\dMx,t)\d^U{x,t)). (3) 

We refer the reader to the classical papers of Gentry, Martin, & Daly [S^] and Harlow & Amsden 
[23] for an interesting discussion on artificial viscosity. Specifically, Gentry, Martin, & Daly [22] 
define the nonlinear viscosity of the type ([s]) to be artificial viscosity, and show that the linear 
viscosity ^ , scaled by the magnitude of local velocity, arises as truncation error (in finite-difference 
approximations). The latter is responsible for stabilizing the transport of sound waves, while ^ 
stabilizes the steepening of sound wavesj^ 

We are primarily concerned with the steepening of sound waves, and shall term artificial viscosity 
of the type ([3]) as classical artificial viscosity. Formally, the use of ^ produces the required amount 
of viscosity near shocks but allows for second-order accuracy in smooth regions. On the other hand, 
the diffusion coefficient \dxu{x,t)\ is precisely the quantity which loses regularity (or smoothness) 
near shock discontinuities. Also, the constant /3 must be larger than one to control numerical 
oscillations behind the shock wave, which in turn overly diffuses the waves and produces incorrect 
wave speeds. 

Alternative procedures have been proposed. For streamline upwind Petrov-Galerkin schemes 
(SUPG), Hughes & Mallet [H] and Shakib, Hughes, & Johan \^S\ use residual-based artificial viscos- 
ity. Guermond & Pasquetti [53] present a similar, entropy-residual-based scheme for use in spectral 
methods. Persson & Peraire [37] develop a method based upon decay of local interpolating poly- 
nomials for discontinuous Galerkin schemes. Later, Barter & Darmofal [281 use a reaction-diffusion 
equation to provide a regularized variant of this approach. 

Our approach is similar to [28j in that it uses a reaction-diffusion equation to calculate a smooth 
distribution of artificial viscosity. Instead of regularizing a DG-based noise-indicator that allows 
for the growth of viscosity near shocks, we regularize the classical artificial viscosity of [21], using 
a gradient based approach for this source term. This approach yields both a discretization- and 
PDE-independent methodology which can be generalized to multiple dimensions by regularizing a 
similar viscosity to that in Lohner, Morgan, & Peraire [25] . 

In 1-D, our approach proves to be a simple way of circumventing the need for characteristic or 
other a priori information of the exact solution to remove oscillations in higher-order schemes. Due 
to the simple and discretization-independent nature of our method, we expect our methodology to 
be useful for a wide range of applications. 

1.3. Outline of the paper 

In Section [2| we introduce the C-method for the compressible Euler equations in one space 
dimension. We show that the C-mcthod is Galilean invariant and that solutions of the C-method 
converge to the unique weak solutions of the Euler equations in the limit of zero mesh size. We also 
show the relative smoothness of our new viscosity coefficient with respect to the classical artificial 
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viscosity of Richtmyer and Von Neumann, and we demonstrate the ability of the C-method to remove 
downstream oscillation in slowly moving shocks. 

In Section |3j we give a brief outline of the numerical schemes whose solutions are used in this 
paper. First, we outline a second-order, continuous Galerkin finite-element method. Second, we 
outline a simple WENO-based finite-volume scheme, only upwinding via the sign of the velocity 
(no Riemann-solvers or characteristic decompositions in primitive variables). The resulting schemes 
applied to the C-method are referred to as FEM-C and WENO-C, respectively. Third, we outline the 
central- finite-difference scheme of Nessyahu and Tadmor (NT), a simple scheme, easily generalizable 
to multi-D [30]. Like our FEM-C scheme, the NT-scheme is at best, second-order, and does not 
require specialized techniques for upwinding. Fourth, we outline a Godunov-type characteristic 
decomposition-based WENO scheme (WENO-G) developed by Rider, Greenough & Kamm [3T] 
which utilizes a variant of a Godunov/Riemann-solver as upwinding, providing a very competitive 
scheme for modeling the collision of very strong shocks. 

In Section |4j we consider the classical shock-tube problem of Sod. With the Sod shock problem, 
we apply our FEM-C scheme and compare with the classical viscosity approach. We then compare 
the FEM-C scheme with the two standalone methods, NT and WENO-G. 

In Section [5j we consider the moderately difficult problem of Osher-Shu, modeling the interac- 
tion of a mild shock with an entropy wave. We compare FEM-C to NT and WENO-G in which 
the differences are more significant than in the Sod-shock comparisons. We show that WENO-C 
compares well with WENO-G; on the other hand, the simple WENO scheme without the C-method 
and without the Gudonov-based characteristic solver also does well in modeling the Osher-Shu test 
case. 

In Section [6j we consider the numerically challenging Woodward-Colella blast wave simulation, 
which models the collision of two strong interacting shock fronts. Though the FEM-C scheme 
performs better than NT, both second-order schemes are somewhat out-performed by the higher- 
order WENO-G method (with characteristic solver). On the other hand, WENO-C compares well 
with WENO-G, having slightly less damped amplitudes with the same shock resolution. 

Finally, in Section [Tj we consider the Leblanc shock-tube, an extremely difficult test case con- 
sisting of a very strong shock. For this problem, devise two strategies to demonstrate the use of 
the C-method. In the first strategy, we use our simplified WENO-C scheme with a right-hand side 
term for the energy equation that relies on a second C-equation which smooths gradients of E/p. 
We obtain an excellent approximation of the notoriously difficult contact discontinuity for internal 
energy, while maintaining an accurate shock speed; simultaneously, we avoid generating large over- 
shoots at the contact discontinuity, which would indeed occur without the use of the C-method. For 
our second strategy, we show that WENO with the Lax-Friedrichs flux can be significantly improved 
with the addition of the C-method. We call this algorithm WENO-LF-C, and show that by using 
just one C-equation (as we have for all of the other test cases), we can sharply resolve the contact 
discontinuity for the internal energy, with accurate wave speed, and without overshoots. 

2. The C-method 

We begin with a description of the I-D compressible Euler equations, written as a 3x3 system of 
conservation laws. We then explain our parabolic regularization, which we call the C-method. 
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2.1. Compressible Euler equations 

The compressible Euler equations set on a 1-D space domain I C M, and a time interval [0, T] 
are written in vector-form as the following coupled system of nonlinear conservation laws: 



u(a;, 0) = Uo(x), x^X,t — 0, 
where the 3-vector u(x,t) and flux function F{u(x,t)) are defined, respectively, as 

and F(u) = 



(4a) 
(4b) 




(E+p) 



and 



Pa{x) 
Uo{x) = I mo{x) 



Eoix) 

denotes the initial data for the problem. The variables p, m, and E denote the density, momentum, 
and energy density of a compressible gas, while p — 'H{p,m,E) denotes the pressure function. It is 
necessary to choose an equation-of-state 'H{p,m, E), and we use the ideal gas law, for which 



p=(7-l) U 



2p) ' 



(5) 



where 7 denotes the adiabatic constant. The equations Q are indeed conservation laws, as they 
represent the conservation of mass, momentum, and energy in the evolution of a compressible gas. 
The velocity field u{x,t) is obtained from momentum and density via the identity 



m 



Inverting the relation ([s]) , we see that the energy density E is a, sum of kinetic and potential energy 
density functions: 



E = 



£?r p 
7-1 

kinetic potential 

The gradient (or Jacobian) of the flux vector F(u) is given by 





(7-3)m^ 

-7^ + (7-1)^ 



7_E 



1 

(3— 7)m 

-(1-7) 



3m/ 
2p2 



with eigenvalues 



\i = u + c, X2 = u, X3 = u 



7-1 

■yni 
P 



(6a) 



where c denotes the sound speed (see, for example, Toro [32,). These eigenvalues determine the 
wave speeds. 
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The behavior of the various wave patterns is greatly influenced by the speed of propagation; as 
such, we define the maximum wave speed to be 



[S'(u)](t) = max max{|Aj(x,<)|} . m 

2—1,2,3 x^X ' — 

as a function of t. 

We are interested in solutions with shock waves and contact discontinuities. The Rankine- 
Hugoniot (R-H) conditions determine the speed s of the moving shock discontinuity, as well as the 
speed of a contact discontinuity. For a shock wave discontinuity, the R-H condition can be stated 

F(U;) - F{\lr) S(U/ - Mr) 

where the subscript I denotes the state to the left of the discontinuity, and the subscript r denotes 
the state to the right of the discontinuity. In general, the following three jump conditions must hold: 

mi-nir ^ s{pi - pr) 
Pi ^ Pi J 

There can be non-uniqueness for weak solutions that have jump discontinuities, unless entropy 



conditions are satisfied (see the discussion in Section 2.9.4 1. So-called viscosity solutions xims are 
known to satisfy the entropy condition (and are hence unique) and are defined as the limit as e — > 
of a sequence of solutions u*^ to the following parabolic equation: 

5tU^ + a,F(u^) =e9,,u% i > 0, (7a) 
u'^=uo, t = Q. (7b) 

In the isentropic setting, for bounded initial data Uq with bounded variation, solutions u*^ converge 
to the unique entropy solution Uyis of Q as e — (see DiPerna [33^ and Lions, Perthame, & 
Souganidis [34]). For non-isentropic dynamics, the same result holds if the initial data has small 
total variation (see Bianchini & Bressan [35]). Moreover, if the initial data Uq is regularized, then 
solutions to ([t]) are smooth in both space and time, and the discontinuity is approximated by a 
smooth function, transitioning from the left-state to the right-state over an interval whose length is 
0(e). 

Some of the classical finite-differencing schemes, such as the Lax-Friedrichs discretization, is 
dissipative to second-order and effectively behaves as a discrete version of Q. The uniform nature 
of such diffusion does not distinguish between discontinuous and smooth fiow regimes, and thus 
adds unnecessary dissipation in regions of the wave profile which do not require any numerical 
stabilization. Such uniform dissipation contributes to a non-phyiscal damping of entropy waves, 
over-diffusion and smearing of contact discontinuities, and changes the wave speeds. Ultimately, 
uniform artificial viscosity is not ideal; rather, artificial viscosity should be added in a localized and 
smooth manner. 
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2.2. Classical artificial viscosity 

The idea of adding localized artificial viscosity to capture discontinuous solution profiles in nu- 
merical simulations dates back to Riclitmyer [T5] , Von Neumann & Riclitmyer , Lax & WendrofF 
[20j . Lapidus |21j and a host of other reseachers. The idea behind classical artificial viscosity is to 



refine the uniform viscosity on the right-hand side of equation ( 7a ) with 



dtu'+d^Fiu')^l3e^d^i\d^u'\d^u^), t>0, (8) 

for a suitably chosen constant /3 > 0, which may depend upon the numerical discretization scheme. 

When the velocity u exhibits a jump discontinuity (i.e., at a shock), the quantity ISajM*^! is 0{j); 
however, away from shocks, where the velocity is smooth, l^^ju'^l remains uniformly bounded in e, 
and in such smooth regions, ([8| adds significantly less viscosity than (7a). On the other hand, as 



we shall demonstrate in Figure [T] the use of l^ajU*^! as a coefficient in the smoothing operator, can 
lead to spurious oscillations in the solution, caused by the lack of regularity in the quantity 

Formally, the use of the localizing coefRcient l^^u'^l corrects for the over-dissipation of the uniform 
viscosity in ([t]) , and a variety of numerical methods have employed some variant of this idea, achieving 
methods that are nominally non-oscillatory near shocks while maintaining second-order accuracy 
away from shocks. However, as we have already noted, the quantity may become highly 

irregular near shock discontinuities, and may thus require setting the constant /? 3> 1 in order to 
stabilize incipient numerical oscillations (see Section |4] for evidence to this observation). While this 
increase in /3 does not effect the asymptotic accuracy of the scheme, it is clearly beneficial to take (3 
as small as possible to preserve the correct amplitude and wave speed. 

The loss of regularity of the coefRcient |i9a;M^| suggests that a smoothed version of IS^u*^! would 
greatly benefit the dynamics. Smoothing l^ajU*^! in space is not sufficient, as we must ensure smooth- 
ness in time as well. Hence, we propose our C-method, which indeed provides a regularized version 
of ^ and allows for the use of much smaller values of (3 (less localized artificial dissipation) , higher 
accuracy, and practical viability. 

2.3. C-method for compressible Euler 



Analogous to ([8|), we control the amount of viscosity in (7a) by the use of a function C{x, t) of 
space and time. This function C{x^t) is a solution to a reaction-diffusion equation, coupling to the 
evolution of u*^. The mechanism of diffusion, smoothing/spreading the sharp peaks localized around 
discontinuities in the velocity competes with the mechanism of reaction, pushing C{x,t) to zero. 
Subsequently, our C-method yields a smooth, yet sharp, distribution of artificial viscosity yielding 
regularized u*^ similar to ([t]) on which we can build high-resolution numerical schemes. 

For fixed Uq we choose /? > 0. Then, for each e > 0, we find 

/ p^{x,t) \ 
u"(a;,t)= nf[x,t) and C"(a;,t) 
V I 

as solutions of the following parabolic system of (viscous) conservation laws: 

dtu' + d^¥{Mi')^d^(pe''C''^d.,M'Y i>0, (9a) 

dtC - S{u')dlC' + ^^C' = S{u'')G{d^u') , i > 0, (9b) 
(u%C^) = (u^,G(9,ii^)), t = 0, (9c) 
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where C^' = + (5 for a fixed positive constant < 6 < Ax, and /3 = ^ maxC^ • forcing to 

X 

equation ^jp) is defined as 

(u^) is defined by (|6]), and Uq denotes a regularization of the initial data which we discuss below. 

max 

We also note that the scaling factor in f3, given by , is included only to make comparisons 

with the classical artificial viscosity approach, but is in no way necessary. 

2.4- Regularization of initial data for use with FEM-C 

Unlike numerical algorithms which advance cell-averaged quantities, the finite-element method 
relies upon polynomial interpolation of nodal values, and requires solutions to be continuous across 
element boundaries in order for the interpolation to converge. As such, the use of discontinuous 
initial data produces Gibbs-type oscillations, at least on very short time intervals. To avoid this 
spurious behavior, it is advantageous to smooth discontinuous initial profiles. 

More specifically, we provide a hyperbolic-tangent smoothing for initial data Uq for our FEM-C 
scheme. Since pointwise evaluation is well-defined for smooth functions, the finite-element discretiza- 
tion scheme can interpolate the regularized data and generate appropriate initial states. 

For an interval [a, 6], we denote the indicator function 

and consider initial conditions with components of the form 

Li 

(u^(x))^ = ^i[,.,,.](x)/;(x), 

where the are pairwise-disjoint (in i) and each /j are smooth. 

We then define the regularized initial condition 



«(a;))^=^lf^.,.,(x)/;(x), 

where 



x + b)\ x + a] 
tanh ^ - tanh ^ 



This regularization procedure achieves approximations of exponential-order away from discon- 
tinuities; near discontinuities, it is a first-order approximation, when measured in the L-'^-norm. 
Specifically, if (uq)* is smooth in a; C I, then the _L^(ci;)-norm of the error 



||(uo)^-«)MIlim= / (uo(x))'-(u^(x))' dx = 0{en (12) 
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for any positive integer p. Alternatively, if Ug is discontinuous somewhere in C X, the L-'^(r2)-norm 
of the error 

||K)-KriU.(o) = 0(e). (13) 

These observations assert that our regularization of the initial data allows for higher-order ap- 
proximation of the initial data and is analogous to the averaging procedure required by Majda & 
Osher \U. 



2.5. A compressive modification of the forcing G in the C -equation 



The function G in (10 1 is chosen in such a manner so that is large where there are sharp 
transitions in the velocity field u'^{x,t). In compressive regions (i.e., where dxU < or dxu'^ < 0), 
sharp transitions over lengths of 0{e) correspond to shocks and artificial viscosity is required so that 

remains smooth. In expansive regions, corresponding to rarefactions, artificial viscosity is not 
generally necessary. 

These observations motivate the following alternative forcing function: 

GcompidxU") = ^ ;t1(_oo,0) {dxu") (14) 

max lOxW^] 

where the indicator function l(_oo,o) introduces viscosity only in regions of compression. 

The ability to use such a switch is heavily dependent on the use of a space-time smoothing. 
Since the velocity in many numerical schemes may become oscillatory near shocks, such a switch 
can become discontinuous between adjacent cells/elements. However, the space-time nature of the 
C-equation resolves this issue, providing a smooth artificial viscosity profile. 

This modified function Gcomp typically increases accuracy in Euler simulations, but can lead to 
a loss of stabilization. For our FEM-C approach, where the stabilizing effects of artificial viscosity 
are necessary to dampen noise, the use of Gcomp is restricted to the problems of Sod and Osher-Shu, 
which contain only moderately strong shocks. 

2.6. Moving to the discrete level 

The use of the C-equation yields smooth solutions u"^ and thus we expect that a variety of 
higher-order discretization techniques, with sufficiently small At and Aa;, could provide accurate, 
non-oscillatory approximations. In our implementation, artificial viscosity spreads discontinuities 
over regions of size 0(/3e). Thus, given a particular initial condition, final time, discretization 
scheme, etc., we choose /3 > such that the scaling e = Ax produces non-oscillatory profiles. 



We also note that the initial condition for C, given in ( 9c I is chosen so to guarantee the coefficients 



of diffusion in (9a) are smooth up to t = 0. Moreover, choosing such initial conditions for C allows 
one to recover the classical artificial viscosity as e — 0. As stated, these initial conditions may 
require a smaller time-step (by a factor of 10) for the first few time-steps. In practice, taking C = 
is an effective simplification to eliminate the need for smaller initial time-steps. Alternatively, we 
can solve an elliptic PDF for C at the initial time and similarly eliminate that concern. 
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2. 7. The C -method under a Galilean-transformation 

We begin our discussion for the case of constant entropy. The Gahlean invariance of the isentropic 
Euler equations results from the advective nature of the PDE. Since we solve a modified equation 
(coupled with the additional C-equation) it is of interest to know to what extent Galilean invariance 
is preserved. For simplicity, we assume that 

p{x,t) = p{x,tf . 

(The choice 7 = 2 corresponds to the shallow water equations, but any other choice of 7 > 1 works 
in the same fashion.) 

Given a fixed w G M we define the change in independent variables 

X = X ~ vt, t = t, 

denoting (f>(x,t) — {x,t) and the analogous change in the dependent variables 

i5{x,t) — p{x vt,t), u{x,t) = u{x -\- vt,t) ~ V. (15) 

A simple calculation yields 

dfp + diipu) = [dtp + dx{pu)] o (f), (16a) 
d^ipu) + diipv? +p)= [dt{pu) + dxipu^)] ocjy + dip-v [dtp + d^ipu)] o <j). (16b) 

We further have that 

p{x,i) = p{x -\- vi,i), (17) 

so that the mass and momentum equations are, in fact, Galilean invariant in the absence of artificial 
viscosity. 

With the C-method employed, ( 16 ) transforms to 

diP + diipu) = [dx{Cdxp)] o (j), (18a) 
diipu) + di{pu^+p) = {dACdM^ - v)]}) o cb, (18b) 

where we let C = e^fiC, and drop the e superscript for notational convenience. 

Examining (9b I, we see that the equation for C is not Galilean invariant, but this is not a physical 
quantity, but can rather be viewed as a parameter to the modified system of conservation laws. As 
such we define the behavior of C under Galilean transformations as follows: 

C(i, t) ~ C[x + wi, t). 

With this definition of C, we find that 

dip + di{pu)= dxiCdip) , (19a) 
diipu) + d.,{pu^ +p) = [dSds,{pu)]} , (19b) 
and hence the C-method for isentropic Euler retains the Galilean invariance. 
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We remark that in the absence of artificial viscosity on the right-hand side of the mass equation, 



the artificial flux term in the momentum equation is modified according to ( 34 ) below. This modi- 
fication ensure Galilean invariance when the mass equation is left unchanged, which is the strategy 
employed for our WENO-C scheme. 

Next, since the Galilean symmetry is for the smooth solutions (for which classical derivatives 
are well-defined), and since smooth velocity fields simply transport the entropy function, it is thus 
a consequence of the transport of entropy, that Galilean invariance holds for the non-isentropic case 
as well. The important of a numerical approximation to capture the Galilean invariant solution is 
fundamental to the initiation of the Kelvin-Helmholtz instability and other basic instabilities present 
in the Euler equation; see Robertson, Kravtsov, Gnedin, and Rudd [33] for a thorough discussion. 
In this connection, we next examine long wavelength instabilities which can arise for very slowly 
moving shock waves. 

2.8. Regularization through the C- equation 

It is of interest to examine the relative smoothness of C to its unsmoothed counterpart \ua:\, 
and to determine the effect of this smoothing relative to the classical artificial viscosity approach. 
In Figure [l] we provide two plots demonstrating the effect of the C-method. In Figure 1 (a) we see 
that the C-equation provides a smoothened viscosity profile compared to the classical approach. 
Alternatively, in Figure 1(b) we plot C using the compression-switch modification Gcomp versus 



using purely the quantity Gcomp (not smoothed by the C-equation) as a viscosity. In both cases 
we see how the C-method provides a far smoother profile with roughly the same magnitude as the 
non-smoothened approach. 




(a) The C-equation with G{ux) (b) The C-equation with Gcomp{ux) 

Figure 1: A comparison of the artificial viscosity profile produced by the C-method and the classical Richtmyer-type 
approach for the Sod shock tube at t = .2. Figures (a) and (b) are with the compression switch off and on, respectively. 
The smooth solid line is the C-method solution, while the oscillatory dashed line is the |na; |-Richtmyer-type viscosity. 

A useful feature of the C-method is the ability to tune parameters in the C-equation to generate 



non-oscillatory behavior. Though we are quite explicit on the form of the C-equation in ( 9b ) , a 



simple modification allows for the diffusion coefficient to be problem dependent, i.e. allowing for a 
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choice of positive constant 7 > and replacing the diffusion term with 



-tS{n^)dlC^ . 

In most of the forthcoming experiments, we fix 7 = 1, but we note that choosing larger 7 can yield 
smoother solution profiles as the profile of C will be less localized. The parameter 7 is a time- 
relaxation parameter, and can be viewed in an analogous fashion to the time-relaxation parameter 
present in Cahn-Hilliard and Ginzburg-Landau theories. For very slow moving shocks, the time- 
relaxation can be adjusted to scale with the shock speed|^ 

We find this to be an effective approach for the flattening procedure discussed in [Hj for removing 
oscillations that form to the left of a slowly right-moving shock. Moreover, Roberts [37] concludes 
that a differentiable form of the numerical flux construction appears necessary to remove downstream 
long-wavelength oscillations caused by slow shock motion. We use the C-method to analyze this. 
Using the slow-shock initial conditions outlined in Quirk [T^], in Figure |2] we show the success 



of the FEM-C (outlined below in Section 3.2 1 in removing these oscillations when choosing 7=1 
(Fig. [2(i^ and 7 = 100 (Fig. [2(b)| ). 




(a) 7 = 1 (b) 7 = 100 

Figure 2: Application of FEM-C to a very slowly moving shock 



2.9. Convergence of the C-method in the limit of zero mesh size 
2.9.1. The isentropic case 

We sketch the proof for the isentropic Euler equations given by 



{pu)t + {pu^ +p)^=0, (20a) 
Pt + (pu)^ = , (20b) 
p{p) = p-^ , (20c) 



^We note that 7 is inversely proportional to the Mach number and its precise functional relation shall be examined 
in future vfork. 
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where 7 > 1. 

To simplify the notation, we set (5 — 1, and set the momentum m'^ = p'^u'^ 



write the C-method version of ( 20 ) as 



Ct - S{u^)Cl, + = S{u^)G{u%) , 



Following ([9|, we 

(21a) 
(21b) 
(21c) 

(21d) 



(m%p^) and flux f(u'=) = ((m^)V/o' + 



or (21i,b) can equivalently in terms of the vector u*^ 
{p'^y , m^) as 

u% + f (u^)^ = e2 [Cu% _ 

where C denotes a diagonal 2x2 matrix with entries C^'^ which is strictly positive-definite. Recall 
that G{u'^) = \u%\/ max satisfies G > 0, and that S{u'^) = maxdu*^ -|-c|, Im*^ — c|), with c denoting 
the sound speed. On any time interval [0,T], the maximum wave speed S{u'^) is uniformly strictly 
positive; thus, as the initial data for C^^q > 0, the maximum principle shows that C'^{x,t) must be 
non-negative. We remark that while the use of C"^'* = C"^ -I- ^ as the coefficient is not required for 
the numerics, as S is taken much smaller than the mesh size Ax, strict positivity of C simplifies the 



proof of regularity of solutions to (21) as well as the convergence argument. 



To avoid issues with spatial boundaries, we shall assume periodic boundary conditions for 
our spatial domain. Note that in this case, the fundamental theorem of calculus shows that 
^ f p{x, t) dx ~ and that mass is conserved. 



2.9.2. The basic energy law 

In order to prove that solutions to (21) converge to solutions of (20), we must establish e- 
independent estimates for solutions of (21). To do so, we multiply equation (21 1) by u"^ , integrate 



over our spatial domain, and make use of the equation (21 d) to find that any weak solution to (21 ) 
must verify the basic energy law 



d 
dt 



lp^{u^rdx+^ 
2 7^1 



dx 



C'^p' [uiydx-e^j 



C^'\pr~\plfdx. (22) 



(The inequality in ( 22 ) is due to the lower semi-continuity of weak convergence and is replaced 



with equality for solutions which are sufficiently regular.) Thus, the total energy of isentropic gas 
dynamics is dissipated according to the right-hand side of ([22]), and for each e > 0, we see that u% 
and p% are square-integrable (in L^) for almost every instant of time, if the density p*^ > A > 0, that 
is, if p'^ avoids vacuum. We shall explain below that this is indeed the case. 



2.9.3. Regularity oj solutions 

The reaction-diffusion equation (21 1) is a uniformly parabolic equation. According to the energy 
law, and as a consequence of Sobolev's theorem, m"^ is a bounded function; furthermore, the right- 
hand side of ( 21 1) is in . It is standard, from the regularity theory of uniformly parabolic equations. 



that then has two spatial (weak) derivatives which are square-integrable. This, in turn, shows 
that for e > 0, solutions possess three spatial (weak) derivative which are square-integrable for 
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Furthermore, by using the symmetrizing matrix 



we can show that t),p'^{-, t)) 



almost every instant of time. This implies that solutions u*^ are classically differentiable in both 
space and time. 

j{pr-') _ 

are, independently of e and t, uniformly bounded in the Sobolev space (consisting of measurable 
functions with two weak derivatives in L^), and thus we may take a pointwise limit of this sequence 
as e — > 0, in the event that the time-interval is sufficiently small as to ensure that a shock has not 
yet formed. Of course, we are interested, in convergence to discontinuous profiles, so we address this 
next. 



2.9.4- Convergence to the entropy solution 

We shall now provide a sketch of the limit as e ^ 0. A function 77 : M is called an entropy 

for ( 20 ) with entropy flux (7 : M'^ — > M if smooth solutions u satisfy the additional conservation law 



r]{u)t +q(u)a; = 0. 



(23) 



In non-conservative form, (pO| and (23) are written as 



ut + V/(u)u^ = 0, V?7(u)ut + Vq(u)u^ = 0, 
from which we obtain the compatibility condition between rj and q, 

Vr;(u)V/(u) = Vg(u). 



(24) 



The pair (77,5) satisfy (23) if and only if condition (24) holds. Moreover, a weak solution to (20) is 
the unique entropy solution if 

77(u)t + <z(u), < . (25) 

For isentropic gas dynamics we can set 

r]{m,p) 



2p 7-1 

which is the total energy, with corresponding entropy flux 

,2 



q{m,p) 



m 
Tp 



7 



7 7 

^P 



V 



m 
P 



We observe that V^rj{m,p) is strongly convex as long as p > 0. 

For the sequence of solution u"^ of (21 1, suppose that as e — 7> 0, u*^ converges boundedly (almost 
everywhere) to a weak solution u of (20). We claim that if (77,(7) satisfy (23), then (251 holds in the 
distributional sense. To see that this is the case, we take the inner-product of V?7(u') with equation 
([2T|), and find that 



77(u^),-|-g(u^), = e2V77(u^)[Cu^L 
Integrating over the spatial domain and then over the time interval [0, T] yields 



ri{\i^{x,T))dx- I 77(u'(x,0))dx 



'■ j[n%fCV\{n^)nldxdt, 
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from which it follows that 







\eul\'^ dx dt < c (26) 



where the constant c depends upon 5, the minimum value of density, and the entropy in the initial 
data. For a smooth, non- negative test function -0 with compact support in the strip I x (0, T), 

jj 'n{u''(t)t+qiu')ct)^dxdt = e J J C(eu')^(/)^ dxdt + J J €^[ulfCV'^rj{u')ul(j)dxdt. 



Thanks to (26 1, the first term on the right-hand side goes to zero like e, while the second term is 
non-negative, since V^?7(u'^) is positive-definite (since 77 is strongly convex) as is C. Thus, as e — > 0, 



we recover the entropy inequality ( 25 ) 



It remains to discuss the assumptions concerning the bounded convergence of u*^ to u, as well as 
the uniform bound from below p"^ > v > 0. The argument relies on finding a priori bounds on the 



amplitudes of solutions to (21 1. If it is the case that uniformly in e > 0, 

|u'| < M and < < , 

then the compensated-compactness approach for isentropic Euler pioneered by DiPerna [33J and 
made much more general by Lions, Pcrthame, & Souganidis [34 provides a subsequence of u*^ 



converging pointwise (almost everywhere) to a solution u of ( 20 ) 



For isentropic gas dynamics, our approximation (21) preserves the invariant quadrants of the 
inviscid dynamics (just as in the case of uniform artificial viscosity) and provides the bound ju*^! < M 
as long as < i> < p"^ for some v. In particular, the Riemann invariants w = u + -^zr^P''^^^ 

z = u — -^zTiP^''^^ satisfy w{x,t) < supit;|t=o and —z{x,t) < sup{—Zt=o) and the intersection of 
these half-planes provides the invariant quadrant (see Chueh, Conley, & SmoUer [3S]), and hence 
the desired bound |u'| < M as long as vacuum is avoided. 

Finally, the fact that we have the lower-bound < < p*^ is an immediate consequence of the 
strong maximum principle. 

2.10. The C -equation as a gradient flow 

Notice that equilibrium solutions to the C-equation are minimizers of the following functional 
(we drop the superscript e): 

£g{C) = J Qc2-G(u.)C+^C2^ dx. 

In the absence of a forcing function G{ux), this reduces to 

UC)^l J (cl + -^C''^ dx. (27) 

The first term is commonly referred to as the Dirichlet energy and its minimizers are harmonic 
functions. The second term can be viewed as a penalization of the Dirichlet energy. In particular, 
because the energy functional is bounded by a constant independent of e > 0, the penalization 
term constrains C to be 0{\/e). Thus, minimizers are trying to be harmonic while minimizing their 
support. 
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The C-equation can be written as a classical gradient flow equation 

where the gradient is computed relative to the inner-product. Thus the heat operator in the 
C-equation, dt — 9^, smooths the forcing in space-time, while the reaction term S{u)/eC minimizes 
the support of the smoothed profile. This is very much related to the theories of Cahn-Hilliard and 
Ginzburg-Landau gradient flows, and we intend to examine this connection in subsequent papers. 



3. Numerical Schemes 

We describe two very different numerical algorithms in the context of our C-method. First, we 
outline a classical continuous finite-clement discretization, yielding FEM-C and FEM-jua;! (based 
on classical artificial viscosity). Second, we discuss a simple WENO-based scheme for compressible 
Euler that upwinds solely based on the sign of the velocity u. To this scheme, we apply a slightly 
modified C-method resulting in our WENO-C algorithm. 

For the purpose of comparison, we also implement two additional numerical methods. The first 
is a second-order central-differencing scheme of Nessayhu-Tadmor (NT), a nice and simple method 
which serves as a base-line for our FEM-C comparisons. The second scheme is a very competitive 
WENO scheme that utilizes a Godunov-based upwinding based upon characteristic decompositions 
(WENO-G). This will serve as a benchmark for our WENO-C scheme. 

3.1. Notation for discrete solutions 

To compute approximations to Q, we subdivide space-time into a collection of spatial nodes 
{xi\ and temporal nodes We denote the computed approximate solution by 

U" W \l{xi,tn), 

noting that for fixed i and n, u" is a 3- vector of solution components, i.e., 

ur = < ■ 

It is important to note that we use the notation u" for both pointwise approximations to u, 
(acquired via FEM-C) and approximations to the cell-average values of u (acquired via WENO-C). 

A subscripted quantity Wi denotes the vector itself and the individual components of the vector. 
We overload this notation so to not cause any confusion between functions defined over a continuum 
versus those defined only at a finite number of points. 

In FEM-C and WENO-C, we discretize ^ (or some slight modification) with e = Aa;, and use 
the above notation for the computed solution. We also denote the approximation to C by C". 
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3.2. FEM-C and FEM-\ux\: A Second-Order Continuous- Galerkin Finite- Element Scheme 

We choose a second-order continuous-Galerkin finite-element scheme to provide a discretization 
of (|9|, subsequently defining our FEM-C scheme. 

We subdivide X with -I- 1 (for N even)-uniformly spaced nodes {xi} separated by a distance 
Aa;. In the FEM community, spatial discretization size is more commonly referred by element-width] 
to maintain consistency with the literature, we refer to the inter-nodal regions as cells. Since we use 
a continuous FEM, the degrees-of-freedom are defined at the cell-edges (as opposed to cell-centers)|^ 
For use in our FEM implementation, it is useful to consider the variational form of ([9]). At the 
continuum level, (u'^,C"^) satisfy 



dt-o." ■ $ - F(u') • 9^$ 



dtC'<i> + S[n')[edxC'dx 



max \dxu'^ 



,2 I 



max C'^ 
I 



dx = , 



e 



dx 



S{u')G{dxu')(t) dx 



(28a) 



(28b) 



for almost every t, for all vector-valued test functions $, and all scalar-valued test functions (p. 

Using the finite-element spatial discretization based on piecewise second-order Lagrange polyno- 
mials, we construct operators -4fem a-nd SfeMj corresponding to the non-time-differentiated terms 



dt 









+ 



= 



(29) 



in (28a I and (28b I, respectively. Using these discrete operators, we write the semi-discrete form of 
(pSafand (|28b|) as 

ApEMiUi, Ci) 
^FEM(Ui, Ci) 

where and Ci represent the nodal values of an approximation to u*^ and C^ for which e = Aa; (see 
Section 2.6 1. For a standard reference on the details of this procedure, see Larsson & Thomee [55] . 

The time-differentiation in ( 29 ) is approximated by a diagonally-implict second-order time- 
stepping procedure; first we predict u""*"^ to and solve the implicit set of equations for C'^'^^ and 

Our fully discrete scheme is given by 



follow by implicitly solving for u" using C 



n+l 



C7 



= C? 



,11+1 



tn + 1 

2 



,cr), 

[^F_EJ\/(u"'''"^, 



\_Afem{'^ 



n+l 



C7 



C7 



^') + s™«,cr)] , 

-')+Afem{K.C^)\. 



(30a) 
(30b) 

(30c) 



For smooth solutions, where artificial viscosity is not necessary, our FEM-C scheme is second- 
order accurate in both space and time when the error is measured in the L^-norm. Moreover, the 
addition the artificial viscosity obtained through the C-method is formally a second-order pertur- 
bation (in Ax) and we have verified this accuracy when /3 > (again, for smooth Uq). For Uq 
containing jump discontinuities, the given scheme is no longer second-order accurate on all of T but 
preserves second-order accuracy in the smooth regions away from discontinuities. 

For the classical artificial viscosity schemes (|8]), the C-equation is no longer used but we require 
a similar step to predict the velocity for use in the diffusion coefficient. This analogous scheme, is 
referred to as the FEM- 1 m,, I scheme. 



^When we compare our FEM-C scheme with other, cell-averaged schemes, we perform an averaging procedure 
based upon averages between nodes. 
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3.3. WENO-C: A Simple WENO scheme using the C -method 

Our WENO-based scheme is motivated by Leonard's finite volume schemes ([3D|, pg. 65). Up- 
winding is performed solely based on the sign of the velocity at cell-edges, and the WENO recon- 
struction procedure is formally fifth-order. 

We divide the interval I into N equally sized cells of width Ax, identifying the N degrees-of- 
freedom as cell-averages over cells centered at the Xi. The cell edges are denoted using the fraction 
index, i.e. 

Xi ~\- Xi-^i 
Xi+1/2 = ^ 

Subsequently, we denote a cell-averaged quantity by Wi and its values at the left and right endpoints 
by Wi_i/2 and Wi+1/2, respectively. 

Given a vector Wi, corresponding to cell-average values, and vectors Zi_i/2, Zi+1/2 corresponding 
to left and right cell-edge values, we define the jth component of vector 

[ WENO (w,,Zj± 1/2)]^. = ^ {w.i + i/2Zj+i/2 - yJj^i/2Zj^i/2) 

where the cell-edge values of Wj^i/2 are calculated using a fifth-order WENO reconstruction, up- 
winding based upon the sign of Zj_^_i^2- 

For the flux in the energy equation, we use 

1 (l + l7) + (l + l^) 



[WEN0£(i?„«,±i/2)]^. = ^(£^, 



J+1/2WJ+1/2 ^ 

Pl^\ ^ J- El 



(l + f^) + (l + li). 

i^.-l/2^.--l/2 (31) 



Using this simplified WENO-based reconstruction, we construct the operators .4weno and 
SwENO where 





Pi 






-^WENO ^ 

















WEN0(p„M,±i/2) 
WENO(m., ^.,±1/2) + dp, ~ ^^"'+V2-gcu.^v. 

WEN0i;(£;„U,±i/2) 



^WENO 



E, 



Ax 



dsCi+i/2 ~ dsC'i^i/2 
Ax 



where for a general quantity Wi, defined at the cell-centers, we denote 



Wi+i/2 = ^ , owi := 

We also use the shorthand notation 

dcUi+i/2 = P Ax"^ max 



-, awi+i/2 ^ 



2Ax 



/2 



max d 



Aa 



Pi+1/2 f^"i+l/25 



(32a) 



(32b) 
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and 



dsC = Ax S{u,) 5Q+1/2. 
Using the above definitions, we define tlie semi-discrete form 





1 


-4wENo(Ui, Ci) 


. ^» . 


^ Ax 


SwENo(Ui,Ci) 



= 



(33) 



and we generate the sequence of iterates u" and C" with a standard fourth-order Runge-Kutta 
time-stepper. 

The resulting discretization outlined above is a slight variation on that outlined in ^ . While the 
amount of artificial viscosity C{x, t) is controlled by only the velocity, we only add artificial viscosity 
to the momentum equation. This change is based upon the fact that WENO already minimizes the 
production of numerical oscillations and the addition of artificial viscosity is primarily intended on 
stabilizing the solution near strong shocks, whereas standalone WENO may lose stability. Without 
dissipation on the right-hand side of the mass equation, it is necessary to modify the artificial 
viscosity on the momentum equation as follows: 



^jid^[Cd^{pu)) ^ e^fid^iCpd^u) . 



(34) 



This modification allows the C-method to maintain a basic energy law (in fact, it is the energy law 



(22 1 with the last term on the right-hand side), and simultaneously permits higher accuracy for our 



WENO-based scheme. 



3-4-. NT: A Second-order Central-Differencing scheme of Nessayhu-Tadmor 

The central-differencing scheme of Nessyahu and Tadmor is an extension of the first-order Lax- 
Fredrichs finite difference scheme in which linear, MUSCL-based reconstructions are used to yield a 
second-order accurate scheme. The resulting scheme is extremely easy to implement (a FORTRAN 
code for 2-D problems is given in the Appendix of |16j ) and does not require the use of Riemann 
solvers or characteristic directions for the purpose of upwinding. The NT scheme allows for various 
choices of limiters to enforce TVD or ENO but the UNO-limiter (see Harten & Osher |3T]) is the 
most successful for our range of experiments. 

Though that NT is easy to implement and is easy generalized to multi-D (yielding the JT- 
scheme [IS]), it merely serves as a base-line comparison for our FEM-C. Both FEM-C and NT are 
second-order, but FEM-C turns out to be far less diffusive by comparison. 

3.5. WENO-G: WENO with Godunov-based upwinding 

In [31] the authors study a fifth-order, WENO-based discretization, upwinding by virtue of a high- 
accuracy Godunov-scheme. Their scheme has the usual trait of WENO, offering minimal diffusion 
near extrema, and has the added stabilization and accuracy of higher-order Godunov solvers. For a 
more in-depth description, see |31j . 



4. Sod shock-tube problem 

For the classic Sod shock-tube problem, we consider the domain T — [0, 1] along with the initial 
conditions 

Po{x) \ / 1 \ / 0.125 \ 

mo(x) = l[o,i)(2:)+ l[i,i](a:), (35) 
E^{x) / V 2.5 / V 0-25 / 
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Figure 3: Comparison of FEM-C and FEM-\ua:\, for the Sod shock-tube experiment with N = 100, T = 0.2. 13 = 0.5 
for both FEM-C and FEM-\u^\. 



imposing natural boundary conditions at x = and x — 1. This standard test problem, first 
considered in Sod [42] . is a preliminary test for the viability of numerical schemes. An exact solution 
is known for this problem and consists of two nonlinear waves (one shock and one rarefaction) along 
with a contact discontinuity. 

In Figure [3(b)] we compare the resuhs of FEM-C and FEM-|ua;| at t = 0.2 using N = 100 cells. 
We note that this comparison uses the standard choice of G in ([9| since we are merely concerned 
with the C-equation performing as a smooth version of classical artificial viscosity schemes. Unlike 
comparisons with the schemes based on cell-averages, we compare the nodal values of FEM-C and 
FEM-|itj;|. In this comparison, we choose /3 — 0.5 for both schemes and see that the accuracy of 
both FEM-C and FEM-|uj:| are quite comparable and each scheme resolves the shock in 3 cells. 
However, we notice noise in FEM-|u2,| near the shock. In Figure [3(b) [ this observation is exemplified 
and we see that FEM-C is relatively non-oscillatory by comparison. 

To limit these oscillations generated by FEM-|ua;|, we increase /3 by a factor of 6 and compare 
the resulting density in Figure [4] In Figure |4(b)| we can see a significant loss in accuracy when 
increasing to /3 = 3. Furthermore, in Figure 4(a) we see FEM-|ua:| requires 6 cells to capture the 
shock. 

In Figure |5] we compare the results of the FEM-C scheme versus NT and WENO-G. Each 
simulation is performed with N — 100 and for the FEM-C scheme we choose (3 ~ 0.4 and now use 
G comp (see Sect ion 2.5|. 

Each scheme produces similar resolution of the shock and contact discontinuity, capturing the 
shock in 3 cells and the contact discontinuity in 6 cells. The NT-scheme produces small, smooth, non- 
physical oscillations as the density transitions from the rarefaction to the lower states, and performs 
the worst at the rarefaction. Both FEM-C and WENO-G are essentially non-oscillatory and despite 
WENO-C performing slightly better at the rarefaction, the results are virtually indistinguishable at 
the shock and contact discontinuity. 
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(a) Complete density profile (b) Zoom-in at shock 



Figure 4; Comparison of FEM-C and FEM-|uj,|, for the Sod shock-tube experiment with N = 100, T = 0.2. /3 = 0.5 
for FEM-C and /3 = 3.0 for FEM-|u^|. 




(a) FEM-C vs. NT (b) FEM-C vs. WENO 

Figure 5; Comparisons of FEM-C against NT and WENO schemes, for the Sod shock-tube experiment with N = 100 
and T = 0.2. 
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5. Osher-Shu shock-tube problem 

For the problem of Osher-Shu, we consider the domain X = [—1,1] along with initial conditions 
Pq{x) \ / 3.857143 \ / 1 + 0.2sin(57rx) \ 

m^{x) = 10.14185 l[_i._o.8)(a^)+ l[_o.8,i] (^): (36) 

En{x) 1 \ 39.1666 / \ 2.5 / 

imposing natural boundary conditions bX x — —\ and x —\ 




(a) FEM-C vs. NT (b) FEM-C vs. WENO 

Figure 6: Comparisons of FEM-C against NT and WENO schemes, for the Osher-Shu shock-tube experiment with 
N = 200 and T = 0.36. 

This moderately difficult test problem, first considered in [3], proves to be more difficult for 
numerical schemes due to the evolution a shock-wave which interacts with an entropy-wave; care is 
required to accurately capture the amplitudes of the post-shock entropy waves. Since the density is 
not monotone, standard flux limiters may unnecessarily apply too much dissipation at local-extrema, 
significantly reducing accuracy. An exact solution for this problem is not available and our 'Exact' 
solution in our plots is generated using the DG-solver furnished in Hesthaven & Warburton [33] with 
3200 cells. 

In Figure [6] we compare the results of FEM-C (we choose (3 — 0.5 and use Gcomp), versus NT 



and WENO-G at t = 0.36. In Figure 6(a) we see that NT diffuses the post-shock amplitudes and 
FEM-C provides far superior results. On the other hand, in Figure [6 (b)| we see that all but one of 
the post-shock amplitudes are slightly better for the WENO-G scheme. This insufficiency of the 
FEM-C scheme is not completely surprising as the FEM-C is only formally second-order versus the 
fifth-order accuracy of the WENO-G scheme. 

Noting this insufficiency of the FEM-C scheme, we compare the WENO-G scheme with WENO-C 



in Figure 7(a) and see the WENO-C scheme is more accurate in resolving the post-shock amplitudes. 
This comes at a price however, as we see WENO-G is more accurate in the N-wave region [—0.6, 0]. 

Furthermore, it is interesting to note that in Figure [7(b) | where we choose /3 = in our simplified 
WENO-scheme, we see that the C-equation is not necessary for Osher-Shu. As we see in Section [6] 
this ceases to be the case as the collision of strong shock waves require stabilization. 



23 




-WENO-C 

- WENO-G 

- 'Exact' 





-WENO-C 
-WENO 
— 'Exact' 



(a) WENO-C vs. WENO-G 



(b) WENO-C vs. WENO 



Figure 7: Comparisons of WENO-C with WENO-G and our WENO scheme with artificial viscosity deactivated, for 
the Osher-Shu shock-tube experiment with N = 200 and T = 0.36. 



6. Woodward-Colella Blast Wave 

The colliding blast wave problem of Woodward-Collella is posed on the domain I — [0, 1] with 
initial conditions 

Po{x) = 1, 
moix) — 0, 

Eo{x) = 250 • l[o.9^i] + 0.25 • l[o.i,o.9) + 2500 • l[o^o.i), 

and reflective boundary conditions at a; = and x = 1. This challenging blast wave problem, 
considered in jl3j tests the ability of a numerical scheme to handle collisions between strong shock 
waves. Any viable scheme generally requires stabilization at these collisions. For the results of a 
wide range of schemes applied to this problem, see [5]. An exact solution for this problem is not 
available and the 'Exact' solution in our plots is generated with a 400-cell PPM solver. 

As is standard in our sequence of experiments, we provide a comparison of FEM-C {/3 = 0.5) 
with NT and WENO-G in Figure[8]at t — 0.038. It is interesting to note, the use of Gcomp is far too 
oscillatory in this difficult test problem; we revert to the standard choice of G. We again see that 
while FEM-C is superior to NT in capturing the amplitude of the two peaks in the density, FEM-C 
is far too diffusive in comparison to WENO-G. 

Despite the relative inefficiency of FEM-C compared to WENO-G, it is interesting to note that 
our FEM-C results (with N = 1200) are better than the artificial viscosity schemes use in Colella & 
Woodward f9|. Our scheme is slightly sharper at the shocks and contact discontinuities and is just 
as accurate in the height of the two peaks. 

Before moving to a comparison of WENO-G and WENO-C, in Figure 9(a)| we see that our 
simplified WENO scheme is highly oscillatory due to the strong shock collision, necessitating the 
use of stabilization. This requirement contrasts to the observations made in Section [5] However, in 
Figure [9 (b)[ we see that the use of a classical artificial viscosity significantly dampens the instability 
but moderate oscillations occur and the C-method provides similar dampening in a smooth way. 
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(a) FEM-C vs. NT 



(b) FEM-C vs. WENO 



Figure 8: Comparisons of FEM-C against NT and WENO schemes, for the Woodward-Colella blast-tube experiment 
with N = 400 and T = 0.038. 





(a) WENO (without stabihzation) 



(b) WENO (with stabihzation) 



Figure 9: WENO with and without stabihzation appUed to the Woodward- Colella blast-tube experiment with N = 400 
and T = 0.038. 
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Finally, in Figure 10 we demonstrate the relative success of WENO-C versus WENO-G. At the 
left peak, WENO-G is more accurate, but at the right peak the reverse situation occurs. Each 
scheme provides very good results, and it is clear that WENO-C is a simple alternative to WENO-G 
which produces similar results for complicated shock interaction. 




Figure 10: Comparison of WENO-C against WENO-G, for the Woodward-Colclla blast-tube experiment with A'' = 400 
and T = 0.038. 





7. Leblanc shock-tube problem 

We conclude our experiments with the Leblanc shock-tube, posed on the domain X = [0, 9], with 
initial conditions 

Po(a;) 

"^o(a;) 1 = 1 I l[o,3)(a;)+ ( |l[3,9](^), (37) 

with natural boundary conditions at x = and a; = 9, and with the adiabatic constant 7 = |. 

Because the initial energy i?o jumps ten orders of magnitude, a very strong shock wave is pro- 
duced, making the Leblanc problem an extraordinarily difficult numerical experiment. First , numer- 
ical methods tend to over-estimate the correct shock speed whenever the shock wave in the pressure 
field is not sharply resolved. Second, numerical approximations tend to produce large overshoots in 
the internal energy 

= P 
' (7-l)p 

at the contact discontinuity. We refer the reader to Liu, Cheng, & Shu [H] and Loubere & Shashkov 
pS] for a discussion of the difficulties in the numerical simulation of the Leblanc problem for a 
variety of numerical schemes. The second-order finite-element basis that we use for our FEM-C 
algorithm is not sufficiently high-order to accurately capture wave speeds in Leblanc, but our fifth- 
order WENO-C scheme is ideally suited for this difficult test case. We shall present two differing 
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strategies for WENO-C, which both capture the correct shock speed and remove overshoots of the 
internal energy. 



7.1. Strategy One: A C equation for the energy density 

As we introduced the C-method in equation Q, artificial viscosity is present on the right-hand 
side of all three conservation laws for momentum, mass, and energy. For the WENO-C algorithm, 
only viscosity in the momentum equation has been used for the Sod, Osher-Shu, and Woodward- 
Colella test cases. Due to the strength of the shock in Leblanc, we now return to using artificial 
viscosity for the energy equation. In our first strategy for this problem, we solve for one additional 
linear reaction-difffusion equation for a new C-coefficient to use on the right-hand side of the energy 
conservation law. 

Specifically, to combat the large overshoot in the internal energy e, we solve a second C-equation 
for the coefficient which we label Cb; the forcing term for the equation uses \dx{E/ p)\/ max 
replacing max \dxu\ which forces the C-equation for the coefficient C„, used for the right-hand 

side of the momentum equation. In particular, since C„ is found using the Gcomp forcing, activated 
only in compressive regions when Ux < 0, for the Ce equation, we activate the right-hand side 
only in expansive regions when Ux > 0. To be precise, this modified WENO-C scheme replaces the 



semi-discrete form ( 33 ) with 



dt 





1 







-4wENo(Ui, Ci) 
SwENo(Ui, Ci) 



(38) 



The resulting fully-discrete scheme solves for u" and 



C" 



where the modified fluxes ^weno and iJwENO are given by: 



A 



WENO 



Pi 

E, 



WEN0(TO„u,±i/2) + ap, 
WEN0£(^,,u,±i/2) 



WEN0(p„u,±i/2) 

3Ct,"i + l/2-t?C„"i 



1/2 



+ 1/2- 



-9cR-Bi-i/2 



(39a) 



WENO 



Pi 

mi 

1 
Ax 



CEi 



-5(u.) 



-^(u,) 



Ce, 



Cu, 
Ce 



C 



comp 



(dui) 



n,d{d{E/ p)i,dui) 



dsCui+^/2 + dsCui_^/^ 
dsCEi^i/2 ~ ^sCe^. 



(39b) 
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The expansive-region forcing for C^; is given by 



G expand i^^ii 



\d{E/p).\ 
ma,x\d{E / p), 



I l[0,oo)(9wi 



(40) 



and we use the shorthand 



dc^Ui+i/2 = Pu Ax^ max dui+1/2 



max Cu 



Pi+1/2 3u.j+l/2 



and 



dcEEi+i/2 = Pe Ax^ max 



+1/2 



C 



Pi+l/2 d{E/p),+ i/2. 



max C Ei ' 



In Figure 11(a) we plot the difference between WENO-C with and without the use of this new 
equation for Ce- For WENO-C with Ce activated, we choose (3^ = 1-0 and Pe = 0.15; with the 
C^-equation deactivated, we use /?„ = 1.0 and ISe = 0. Observe that activating the C^j-equation 
removes the large overshoot at the contact discontinuity. Furthermore, examining the location of 
the shock, we see that the use of the C£;-equation produces more accurate approximations of the 
shock speed. 




(a) With and without Ce, N = 1440 (b) Successive refinements, N = 360, 720, 1440 

Figure 11: Internal energy plots for WENO-C for the Loblanc shock-tube experiment at T = 6. 



In Figure |ll(b)] we show the results of WENO-C at iV = 360, 720, 1440. In this plot, we see very 
little overshoot at each level of refinement and this small overshoot does not grow with refinement. 

7.2. Strategy Two: a new type of viscosity for the energy density 

Our second strategy for the Leblanc problem may be viewed as being motivated by the energy 
dissipation rate of real fluids, and adheres to our framework of only solving one C-equation, forced 
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by the normalized modulus of the gradient of velocity. The idea is easy to explain, and we begin by 
writing the equations for momentum and mass (we drop the superscript e): 



Ct - S{u)C,^, 



Pt + {pu)x = 
Siu) 



-C = Siu)G{ux) . 



(41a) 
(41b) 
(41c) 

(41d) 



By multiplying the momentum equation by the velocity u, integrating over the spatial domain, and 
using the conservation of mass equation, we find the basic energy law: 



1 2 , ^ 1 

-pu ax H 

2 7- 



p dx 



= —e (3 Cpu^ dx . 



(42) 



Note, that when e = 0, the variable E is exactly the energy density; that is, when e — 0, E = 



Ipu" 



T^hi- Thus, we formulate a right-hand side term for the energy equation to ensure the E 



continues to represent the energy density for e > 0. To do, we choose a right-hand side which will 



provide the same energy law as (42). We introduce the following equation: 

{uE- 



Et 



Up)x 



-e^jiCpul 



(43) 



The fundamental theorem of calculus shows that integration of ( 43 1 provides the same basic energy 



law as (42 ). Hence, our second strategy employs the equation (41 ) together with (43 ). The interesting 



feature of the new right-hand side of the energy equation is its nonlinear structure, quadratic in 
velocity gradients. This energy loss compensates for entropy production, and can become anti- 
diffusive near contact discontinuities. As such, we shall discretize this set of equations using the 
very stable Lax-Friedrichs flux. We remark that the term e^pCpu^ is analogous to the viscous 
dissipation term of the Navier-Stokes- Fourier system and can be found as a truncation error in }22| . 

As we noted above, to the best of our knowledge, the most commonly used numerical schemes 
applied to Leblanc tend to exhibit a significant overshoot in the internal energy e at the contact 
discontinuity. Furthermore, on coarse meshes (< 2000 cells), the speed of the shock tends to be 
inaccurate. Indeed, this is the case for arguably the most widely used WENO implementation, 
designated WENO-LF-5-RK-4 by Jiang & Shu [5]. This scheme, which we call WENO-LF, uses a 
Lax-Friedrichs flux-splitting with a 5th-order WENO reconstruction in space and 4th-order Runge- 
Kutta in time. 

If we examine the contact discontinuity at x ~ 6.8 in Figure 12(a) 
we see that WENO-LF exhibits relative overshoots of 12. 



at resolutions N = 360, 720, 1440 



%, 11.8% and 11.4% respectively This 
slow decay of the overshoot suggests that WENO-LF suffers from the Gibbs-phenomenon, despite 
it's attempt to quell oscillatory behavior. Examining the shock at x « 8 we see that the computed 
shock speeds are inaccurate. 

To address the loss of accuracy exhibited by WENO-LF, we propose the use of the C-equation 
along with a nonlinear viscosity on the energy equation. Since WENO-LF has an intrinsic artificial 
viscosity (by virtue of the Lax-Friedrichs splitting) on the right-hand side of the momentum equation, 
we find that we do not need to explicitly use our artificial viscosity for the momentum (even though 
this mathematically motivated our nonlinear viscosity for the energy equation). As such, we require 
a single C-equation which is forced by Gcomp{ux)- 
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Keeping consistent with the semi-discrete formulation, we write the WENO-LF-C scheme 





1 







"4wENO-LF(Ui) + 'H(Ui, d 
SwENo(Ui, Ci) 



(44) 



where ;Sweno is given by ( 32b ) and .4weno-lf corresponds to the choice of the WENO flux described 
in P (i.e. if H = then (44) is the same as WENO-LF). The term H is a discrete approximation 
of /3e^C"^'''/9'|3a;M'p. The operator H is defined as 



H(u„a) 







/3 Ax^ max 1 1 / 2 1 j^^^ 1 



In Figure |12(b)| we demonstrate the benefit of WENO-LF-C with /3 = 5.0, again at successive 
refinements o^ N — 360, 720, 1440. The overshoot at the contact discontinuity is relatively non- 
existent while the shock speeds are far more accurate and appear to converge to the correct speed 
at a faster rate. 



- WENO-LF 360 ceils 

- WENO-LF 720 cells 
-WENO-LF 1440 cells 



- WENO-LF-C 360 cells 

- WENO-LF-C 720 cells 
-WENO-LF-C 1440 cells 



(a) WENO-LF 



(b) WENO-LF-C 



Figure 12: Internal energy plots for the Leblanc shock-tube experiment at T = 6 using WENO-LF with and without 
the C-equation. 



8. Concluding Remarks 

We have presented a localized space-time smooth artificial viscosity algorithm, the C-method, 
and have demonstrated its efficacy on a variety of classical one-dimensional shock-tube problems. As 
compared to more established procedures, the C-method has been shown to be highly competitive 
with regards to accuracy and stability, while being relatively easy to implement. Because of its 
simplicity, the C-method can readily be extended to multiple space-dimensions and/or utilized in 
reactive-fiow simulations. Of value to reactive fiows is the localized smooth diffusion provided by the 
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C-method; specifically, the function C can be used to actively influence various mixing-ratc-limitcd 
reactions occurring near sharp boundaries. 

In the future, the gradient-based source term used in the current implementation of the C-method 
may be combined with a noise-indicator that turns off the current gradient-based source term when 
it is not needed. Such noise-indicators require a very high-order scheme compatible with DG or 
llth-order WENO to name just two examples. By projecting the solution onto a suitable basis, the 
noise-indicator would activate when small-scale coefficients of this basis do not have sufficient decay; 
in turn, an indicator function, localized about the region of noise, would activate and force the C 
equation. This approach is taken in [28], but without any gradient-based forcing functions like our 
function G or Gcomp- 

For example, after the rapid initial growth of the internal energy field in the Leblanc shock-tube 
problem, this field is essentially representative of the advection of a square-wave. Thus, after initial 
growth, the gradient-based source term in the C equation for energy could be deactivated leading 
to less diffusion in the downstream contact discontinuity; simultaneously, the noise-indicator would 
activate if small-scale instabilities were to set in. 

But, for more general problems, the impact of the activation/deactivation of the source term in 
the C-method on numerical accuracy is not entirely obvious and is left for future research. 
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